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ABSTRACT 

The Hi 21 cm transition line is expected to be an important probe into the cosmic dark ages and epoch of 
reionization. Foreground source removal is one of the principal challenges for the detection of this signal. This 
paper investigates the extragalactic point source contamination and how accurately bright sources (> 1 Jy) 
must be removed in order to detect 21 cm emission with upcoming radio telescopes such as the Murchison 
Widefield Array (MWA). We consider the residual contamination in 21 cm maps and power spectra due to 
position errors in the sky-model for bright sources, as well as frequency independent calibration errors. We 
find that a source position accuracy of 0.1 arcsec will suffice for detection of the Hi power spectrum. For 
calibration errors, 0.05 % accuracy in antenna gain amplitude is required in order to detect the cosmic signal. 
Both sources of subtraction error produce residuals that are localized to small angular scales, A:j^ > 0.05 Mpc"' , 
in the two-dimensional power spectrum. 

Subject headings: Cosmology: Early Universe, Galaxies: Intergalactic Medium, Radio Lines: General, Tech- 
niques: Interferometric, Methods: Data Analysis 



1. EVITRODUCTION 

The cosmological epoch of reionization (EoR) is a key 
milestone in the history of structure formation, marking the 
transition from a fully neutral to a highly ionized intergalactic 
medium (IGM) due to the ultra-violet and X-ray radiation of 
early stars, galaxies, and black holes. Recent observations of 
the Gunn-Peterson effect, i.e., Lya absorption by the neutral 
IGM, toward the most distant quasars (z ^ 6), and the large 
scale polarization of the CMB, corresponding to Thompson 
scattering during reionization, have set the first constraints 
on the reionization process. These results suggest significant 
variance in both space and time, starting perhaps as far back 
as z -- 1 1 (Komatsu et al. (2010); WMAP 7 year data) and ex- 
tending to z ^ 6 (Fan et al. 2006). Previous WMAP five year 
data indicates the 5 a detection of the E-mode of polariza- 
tion which rules out any instantaneous reionization at z ^ 6 at 
3.5 (T level. In case of the Gunn-Peterson effect, the IGM be- 
comes optically thick to Lya absorption for a neutral fraction 
as small as ^ lO""'. In order to overcome these limitations, it 
has been widely recognized that mapping the red-shifted Hi 
21 cm line has great potential for direct studies of the neutral 
IGM during reionization (Barkana & Loeb 2001; Fan et al. 
2006; Furlanetto et al. 2006; Morales & Wyithe 2009). 

There are number of upcoming low-frequency arrays with 
key science goals to detect the Hi 21 cm signal from the 
EoR. This includes the Murchison Widefield Array [MWA] 
(Mitchell et al. 2008; Lonsdale et al. 2009; Bowman et al. 
2006), Precision Array to Probe Epoch of Reionization [PA- 
PER] (Backer et al. 2007; Pai'sons et al. 2009), Low Fre- 
quency Array [LOFAR] (Harker et al. 2010; JeUc et al. 2008; 
Labropoulos et al. 2009) and Giant Meterwave Radio Tele- 
scope [GMRT] (Pen et al. 2009). One of the major challenges 
for all of these upcoming arrays will be the removal of the 
continuum foreground sources in order to detect the faint Hi 
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signal from the EoR. 

A variety of continuum foregrounds complicate redshifted 
21 cm measurements of the EoR (Shaver et al. 1999). Diffuse 
Galactic synchrotron emission dominates the low-frequency 
radio sky and is approximately four orders of magnitude 
brighter than the ~ 10 mK 21 cm signal at the frequency rele- 
vant to reionization (i/ w 150 MHz). In addition. Galactic and 
extragalactic free-free emission contribute additional flux to 
the diffuse foreground. Radio point sources from AGN, radio 
galaxies, and local Galactic sources are numerous and partic- 
ularly challenging. The brightest of these sources have fluxes 
well above 5 > 1 Jy and are seven or eight orders of magni- 
tude above the EoR signal in low-frequency radio maps. The 
distribution of point sources also extends to very faint levels 
such that the brightness temperature due to confused sources 
in upcoming arrays will be ^ 10 K, or three orders of magni- 
tude brighter the than the 21 cm background. 

In this paper we discuss how the radio interferometric imag- 
ing techniques are going to affect the foreground source mod- 
eling and subsequent removal from the data-set in order to 
search for the EoR signal. Recently, there has been exten- 
sive research on foreground source modeling at these low fre- 
quencies (Di Matteo et al. 2002; Jelic et al. 2008; Thomas 
et al. 2009). Similar effort has also been made in exploring 
different techniques to remove the foregrounds from the EoR 
data-set by Morales et al. (2006a,b); Bowman et al. (2009); 
Liu et al. (2009); Parsons et al. (2010); Gleser et al. (2008); 
Harker et al. (2009b, a). Since attempting to observe a sig- 
nal below the confusion limit of foreground sources is a novel 
aspect of 21 cm experiments, most of these works primarily 
focus on the removal of faint and confused sources that fall 
below a specified cutoff flux limit. Sea (~ 1 Jy)- They do not 
consider the foreground sources brighter than Scut and how 
accurately they need to be removed. Indeed, most of these 
analysis implicitly assume that the bright foreground sources 
above Scm have been removed perfectly. But in reality im- 
perfect instrument calibration or any errors in the subtracted 
foreground model will introduce artifacts and leave residual 
contamination in the data after bright source removal, even by 
traditional techniques such as "peeling". These residuals may 



interfere with either the subsequent faint source subtraction or 
the ultimate detection and characterization of the redshifted 
21 cm signal. 

In Datta et al. (2009) we dealt with the bright point sources 
above Scur and the limitations that will be caused due to imper- 
fect removal of such sources in the image plane. In this paper 
we extend the initial analysis in order to estimate the resid- 
ual contamination in the power spectral domain of improper 
bright source subtraction. The objective of this paper is to 
demonstrate how the accuracy in the foreground removal af- 
fects the detection of Hi 21 cm power spectra with the MWA. 
In Section 2, we discuss our choice of sky model and out- 
line the simulation parameters, including the array specifi- 
cations and data reduction procedure, and describe the two 
categories of corruption terms that we will consider: source 
position errors and residual calibration errors. The results 
obtained for the residual angular power spectrum, spheri- 
cally averaged three-dimensional power spectrum, and two- 
dimensional power spectrum are presented in Section 3. Fi- 
nally, in the last section we summarize the implications of the 
results from our simulations. 

2. THE SIMULATION 

2.1. Sky Model 

Our main aim is to explore the level of accuracy needed 
in instrument calibration and foreground modeling in order 
to ensure the residual errors from bright foreground source 
removal do not obscure the detection of the signal from cos- 
mic reionization. With this goal in mind, we use a simple sky 
model for our simulations that only includes bright radio point 
sources. No diffuse emission from the Galaxy is included as a 
part of the sky model; and the 21 cm signal and thermal noise 
are also omitted. Our sky model is derived from the logA^ 
- log S distribution of sources and is termed the "Global Sky 
Model" (GSM) from now onwards. Since the GSM only in- 
cludes sources above 1 Jy, we follow the source counts from 
the 6C survey at 151 MHz (Hales et al. 1988): 



TABLE 1 
Array Specifications 
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For a field-of-view of 15° the total number of sources above 
1 Jy is ^ 170, following the above power-law distribution. 
The entire flux range, between 1 Jy and 10^ Jy, has been sub- 
divided into several bins (in logarithmic scale) and populated 
with the number of sources that corresponds to the flux range 
of each bin (according to Equation 1). Inside each bin, we 
have assigned each source a flux density following a normal 
distribution. The strongest source in our GSM is ^ 200 Jy. 
The observed distribution of radio sources shows evidence 
for only very weak angular clustering and the brightest ex- 
tragalactic sources in the sky are not clustered at all (Blake & 
Wall 2002). Therefore, in order to assign a position to each 
of these sources within the field-of-view, a uniform random 
number generator has been used which predicted the offset 
from the field center for respective sources. In the GSM all 
the foreground sources are flat spectrum, i.e. with zero spec- 
tral index (a = 0). 

Figure 1 shows a simulated image of our bright source 
sky model that has been produced using the procedure de- 
scribed below. A wide-field variant of the well-known Clark- 
CLEAN algorithm that utilizes a w-projection algorithm for 
3-dimensional imaging was applied to the simulated image. 
The apparent angular size of the sources in Figure 1 reflects 
the size of the synthesized beam. The input sources in the 



Parameters 


Values 


No. of Tiles 


512 


Central Frequency 


158 MHz (z ~ 8) 


Field ot View 


~ 15" at 158 MHz. (oc A) 


Synthesized beam 


~ 4.5' at 158 MHz. (oc A) 


Effective Area per Tile 


~ 17 m^ 


Maximum Baseline 


~ 1.5 km 


Total Bandwidth 


32 MHz 


Tsy, 


~ 250 K 


Channel Width 


~ 32 kHz 


Thermal Noise 


~ 7.55 mK 
(5000 hours & 2.5 MHz) 



Note. — Array parameters have been influenced 
by the MWA specifications as mentioned in Mitchell 
et al. (2008) and Bowman et al. (2009). Original MWA 
Field-of-view is ~ 25° at 150 MHz. 

model are treated as ideal point sources. This input sky model 
is used for all the simulations presented in this paper. 
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Fig. 1 . — Simulation of the sky model centered on RA=4h and DEC=-26° 
as would be observed by the MWA. Clark-CLEAN has been applied to this 
image using w-projection (256 planes) and natural weighting (Datta et al. 
2009). 



2.2. Array Specifications 

Table 1 outlines the instrumental parameters that we have 
assumed for this analysis. Most of these parameters reflect the 
current specifications for the MWA, but we note that the array 
is presently under development and some properties may be 
subject to change. In addition, we have intentionally reduced 
the simulated field of view compared to the actual MWA in 
order to reduce the computational overhead of the simulation. 
Figure 2 shows the array layout for the 5 12 element array with 
maximum baseline of 1 .5 km. 

For the purposes of modeling earth rotation synthesis in the 
instrumental response, the center of the target field is chosen 
such that it coincides with one of the cold spots in the fore- 
ground Galactic synchrotron emission visible from the south- 
ern hemisphere location of the MWA. The exact field cen- 
ter used for the GSM is 4 hours in Right Ascension and -26 
degree in Declination. Most of the upcoming low-frequency 
telescopes, including the MWA, will only be able to observe 
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Fig. 2. — Array layout for the 512 elements with maximum baseline of 
1.5 km. 

a field around its transit. We have used 6 hours of integrations 
for all the simulations, assuming that the telescopes will ob- 
serve a field between ± 3 hours in Hour Angle from transit. 

2.3. Data Reduction Procedure 

The 15° field-of-view will include ^170 bright sources 
(> 1 Jy). The individual flux densities of these foregrounds 
are ^ 10^-10^ times higher than the signal from cosmic 
reionization that these instruments are aiming to detect. So 
the challenge lies in calibration and subsequent removal of 
such bright sources from the raw data-sets. The data rate of 
19 GB/s (Mitchell et al. 2008) will not allow the MWA to 
store the raw visibilities produced by the correlator. Hence 
real-time calibration and imaging needs to be done in order to 
reduce the data volume and store the final product in the form 
of image cubes (Mitchell et al. 2008). The critical steps in- 
clude removal of the bright sources above the Sai level from 
the data-sets in these iterative rounds of real-time calibration 
and imaging procedure. As a result the residual image-cubes 
will not be dominated by these bright sources and the rest of 
the foregrounds can be removed in the image domain. 

However, the accuracy of the foreground source removal 
strategies are strongly dependent on the data reduction pro- 
cedure. The likely data reduction procedure which will be 
followed by the upcoming telescopes can be broadly outlined 
as : 

• The raw data-sets from the correlator will go through 
real-time calibration and subsequent removal of the 
bright sources based on some Global Sky Model 
(GSM), down to Scut level, in the UV domain. 

• The residual data-sets will be imaged and stored as a 
cube for the future processing and removal of sources 
which are below Scut- 

The simulated data reduction pathway that we follow in this 
paper is: 



(i) First, the observed visibilities 
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(y,.f' (u,i') = 



(u, i^)) are simulated for a 6 hour observation 
(±3 hours in Hour-angle) using the GSM and the array 
configuration from Section 2.2. In the above notation, 
n^v = (m,v) denotes the Fourier conjugate of the sky 
coordinate (9x,0y) and v is the frequency of observations. 

(ii) Next, we generate the foreground model (y/"°"'(u,i^)) 
that will be subtracted from the observation. In this case, the 



model is corrupted to either simulate errors in the assumed po- 
sitions of the sources or to simulate calibration errors. For the 
source position errors, the model visibilities are given simply 
by: 

jGSMjiiiperfec 



v,.r'(u,^): 



y. 



'(u,iy), 



(2) 



where the position of each source has been slightly moved 
from its original location by a distance drawn from a Gaus- 
sian distribution with standard deviation ag. We assume that 
the source position errors are constant throughout the entire 
duration of the observation, as would be the case for a fore- 
ground model constructed from either an outside catalog or 
from the data itself at the conclusion of the observation. This 
is an idealization that may be broken in practice if sources are 
"peeled" in real-time. 

For the residual calibration errors, the model visibilities are 
given by: 



v!j'"'{u,i^)=g!(t)g*itw; 



GSM,„ 



'(U,^), 



(3) 



where gi(t) w (1 +ai)e"^' are the antenna-dependent complex 
gains. The parameters a, and 0, denote small amplitude and 
phase deviations, respectively, and are each drawn from there 
own Gaussian distribution with standard deviation da or a^ 
(Datta et al. 2009). 

The MWA will produce calibration solutions in real-time 
with an ^ 8 second cadence. This rapid pace is planned in or- 
der to simultaneously calibrate both the instrument hardware 
properties and the ionospheric phase screen. It is not known, 
yet, if in practice the residual calibration errors from the real- 
time processing will be largely independent or highly corre- 
lated between individual 8-second solutions. This is an im- 
portant experimental property to consider in our simulation, 
because the degree of correlation greatly affects the level of 
accuracy needed in individual calibration solutions. If the in- 
dividual errors are largely independent, then each 8-second 
sample can be modeled as coming from a Gaussian distribu- 
tion and the accuracy tolerance will be relatively loose since 
many samples will be available and tend to average toward 
zero. Such a situation would be the best-case scenario. On the 
other hand, if the calibration errors are highly correlated, then 
each calibration solution must meet a much more stringent ac- 
curacy level to achieve the same residual contamination at the 
end of the full integration. 

For our simulation procedure, we assume a relatively con- 
servative scenario that the residual errors in a given antenna's 
8-second calibration solutions are perfectly correlated for the 
duration of one 6-hour observing night, but perfectly uncorre- 
cted between successive observing nights. We further assume 
that the residual errors between antennas are perfectly uncor- 
rected at all times. This choice is somewhat arbitrary given 
the current level of knowledge, but we believe it is a plausible 
fiducial case since both the overall ionospheric properties and 
the ambient conditions may change significantly from day-to- 
day. Hence, in our simulation, Ca and aij, are used to draw a 
calibration error value (a,- and 0,) from a Gaussian distribu- 
tion only once per antenna per night and that specific error is 
applied to all the simulated 8-second solutions for the given 
antenna throughout the 6-hour period of rotation synthesis. 
When the next night's observing block commences, a new er- 
ror is drawn from the distribution for each antenna, and so 
on. 

(iii) Now we are ready to calculate the residual visibili- 
ties by subtracting the foreground model produced in step (ii) 
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Fig . 3 . — (a) The spectral profile along two lines of sight in the the final residual image cube following the IMLIN step. In this case, the two pixels were chosen 
to be next to the positions of sources in the input model, (b) Same as panel (a), but here, the two pixels were chosen to be far from any sources in the input model. 
Synthesized beam area of 4.5 arcmin X 4.5 arcmin is used to convert the flux densities to surface brightness. 

affects the end result of the bright source removal, so we ap- 
proximate the faint source polynomial fitting here by applying 
a Fourier transform to the UV map generated from the resid- 
ual visibilities in order to produce a residual dirty image cube, 
V^ In this dirty image cube, we fit a third order polynomial 
in frequency along each line of sight and subtract it. Thus, we 

obtain the final residual image, I;mun(^^ ^)- ^^ refer to this 
step as "IMLIN" for the remainder of this paper because it 
was implemented with the IMLIN algorithm (Cornwell et al. 
1992). To illustrate this final result. Figure 3 shows residual 
spectral profiles along four lines of sight in the dirty image 
cube after polynomial fitting and subtraction. Example im- 
ages of the final residual contamination after IMLIN are also 
shown in Datta et al. (2009, panels (b) of their Figures 5, 7, 
and 9). 

Using higher order polynomials in the IMLIN step re- 
moves structures at increasingly smaller scales (McQuinn 
et al. 2006). This improves the foreground cleaning, but since 
the 21 cm reionization signal has significant structures on 
scales that correspond to ^2.5 MHz, or approximately 10% 
of the bandwidth over which the polynomial is fit, it also has 
the potential to remove much of the 21 cm signal. We have re- 
stricted our attention to a third order polynomial in this work 
because it is the lowest-order polynomial likely to be suffi- 
cient for removing the faint continuum sources given their 
power-law spectral shapes. 

We also explored using the UVLIN algorithm (Cornwell 
et al. 1992), which fits and subtracts polynomials in the UV 
domain instead of the image domain, eliminating the need to 
convert our residual data sets into image cubes. However, 
UVLIN works perfectly only within a small field-of-view, de- 
pending on the channel width in frequency (Cornwell et al. 
1992), and was found to be inadequate for our purposes. 

(v) The final step in our procedure is to calculate power 
spectra from the residual image cubes and compare these 
residual foreground power spectra to the theoretically pre- 
dicted 21 cm power spectrum and expected thermal noise 
power spectrum for the MWA. We calculate three forms of 
the residual power spectra from our final data cubes: the de- 
rived angular power spectrum Q for a narrow frequency chan- 
nel, the spherically averaged three-dimensional power spec- 
trum P(k) from the entire data cube, and the two-dimensional 



Fig. 4. — ID spherically averaged power spectrum of the input Global Sky 
Model showing the total power of the bright sources in the sky model before 
any foreground removal has been applied. The thermal noise uncertainty for 
a 300 hour observation by the MWA is also shown, along with the Hi 21 cm 
signal power spectrum for a fully neutral IGM (jthi = 1) at z = 8 (Furlanetto 
et al. 2006). 
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from the simulated observation of step (i) according to: 

y,f(u,i^) = y,f(u,i.)-v,7'"'(u,i.). 

In the sections below, we refer to this step as "UVSUB" since 
it was implemented using the UVSUB algorithm (Cornwell 
et al. 1992). For the residual calibration errors, we can reduce 
Equation 4 by substituting in with Equation 3 and simplifying 
to obtain: 
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(iv) At this point, we have completed the subtraction of 
bright sources from our simulated observation, leaving only 
residual contamination due to the differences between our 
simulated observation and the corrupted model. Example im- 
ages of the residual contamination at this stage are shown in 
Datta et al. (2009, panels (a) of their Figures 5, 7, and 9). 
In practice, this bright source-subtracted data cube will be 
the starting point for the second stage of redshifted 21 cm 
foreground subtraction that aims to remove faint and con- 
fused sources by fitting and subtracting a low-order polyno- 
mial along the frequency axis for each line of sight in the 
data cube. We want to understand how this additional process 



power spectrum P(k±,k\0 found by averaging over transverse 
modes in the full three-dimensional power spectrum. Each of 
these cases is discussed in more detail in Section 3. As a ref- 
erence, we show in Figure 4 the spherically averaged power 
spectrum for our input GSM before any source removal has 
been applied. 

In order to simulate the observed visibilities (V,y'"), we have 
used the simulator tool in the CAS A software ^. We have 
also used CASA to perform the imaging and the subsequent 
IMLIN step. The rest of the operations are performed using 
separately written python ^ scripts. 

3. RESULTS 

We begin our discussion of the results of the residual power 
spectrum determination by reviewing our initial findings from 
Datta et al. (2009). In that work, we explored the source po- 
sition and calibration accuracy needed to allow direct imag- 
ing of Stromgren spheres with very deep integrations by the 
MWA. Our simulations demonstrated that knowledge of the 
true positions of the bright foreground sources in an MWA 
target field is required to within ere =0.1 arcsec, assuming 
Gaussian errors, in order for the residual contamination fol- 
lowing subtraction to be below the 21 cm signal from Strom- 
gren spheres in image maps that could be acquired by the 
MWA with 5000 hours of integration. Similarly, in Datta 
et al. (2009) we found that, for the case of calibration er- 
rors corrupting the measurements under the same conserva- 
tive assumptions outlined in step (ii), a calibration accuracy of 
aa = 0.2% in gain amplitude (or a^, = 0.2 degree in phase) is 
needed for the residual contamination to be below the thermal 
noise in a part of the image map far from any bright sources 
for a long integration by the MWA. 

Here, we focus our attention on the residual contamina- 
tion that can be tolerated in measurements of the power spec- 
trum of a target field observed by the MWA. One of the pri- 
mary motivations for seeking to first detect and characterize 
the 21 cm power spectrum, rather than immediately attempt 
to image the background, is that the MWA will only require 
^ 300 hours of observing to have sufficient sensitivity to de- 
tect the spherically-averaged 21 cm power spectrum at z « 8, 
assuming the IGM is not fully ionized at that time. Because 
the power spectrum measurement differs significantly from 
a direct imaging observation, there are several key questions 
that we seek to address: 1) are the tolerances on the source 
position and calibration errors greater (or lesser) than in the 
direct imaging case, 2) is a particular region of the power 
spectrum likely to be more affected by residual contamina- 
tion than another, and 3) can we hope to build a library of 
template models for foreground contamination that could be 
used to marginalize out some of the contamination during the 
analysis of the power spectrum? 

We will address these questions for each of the three classes 
of power spectra listed in step (v) (section 2.3) that the MWA 
is likely to study: angular, spherically averaged, and two- 
dimensional. For each class of power spectrum, we will 
present residual power spectra for both corruption models: 
source position errors, and calibration errors. And for each 
corruption model, we will use three fiducial levels of error in 
our investigation: for source position errors, our fiducial cases 
are a-g = {0.01, 0.1, and 1 arcsec}, while for the calibration er- 
rors our fiducial levels are cTa = {0.01, 0.1, and 1%} in gain 

' http://casa.nrao.edu/ 
' http://www.python.org/ 



amplitudes, which also translates to a-p = {0.01, 0.1, and 1°} 
in phase (Datta et al. 2009). 

3.1. Angular Power Spectrum 

White et al. (1999) describe the technique to derive the an- 
gular power spectrum from radio interferometric data. Using 
the flat field approximation (Datta et al. 2007): 



Q = 



E2 



,W(u)\V(u)\^ 



E2.\u\=iW(u) 



(6) 



where |u| = vt^+v^ and i ~ 27r|u| under flat-field approxi- 
mation. Here, V(u, v) is the un-weighted visibilities from the 
residual images and W(u) is the number of visibilities enter- 
ing each u cell. 

In Figures 5 and 6, we have plotted i{l+ l)Q/(27r) calcu- 
lated for one frequency bin of width 125 KHz from our simu- 
lated residual image maps. Figure 5 shows the angular power 
spectrum resulting from using the foreground model that is 
corrupted by source position errors. Figure 6 illustrates the 
same result for the case of residual calibration errors. The 
top panels (a) of both figures show the angular power spectra 
derived from the UVSUB residuals of step (iii) in our analy- 
sis procedure. The bottom panels (b) show the angular power 
spectra from the final IMLIN residuals following step (iv). 
The vertical lines in panel (b) of both figures corresponds to 
(. ^ 250. The IMLIN step, which involves fitting a third order 
polynomial over a total bandwidth of 32 MHz, is expected to 
remove most of the significant structures for scales larger than 
this (corresponding to £ < 250). All of the plots have been re- 
stricted to i* < 5000 to match the size of the MWA synthesized 
beam (4.5 arcmin). 

The total thermal noise power is much stronger than the 
angularpowerspectraof the Hl21 cm signal. Hence, we have 
assumed that the final power spectrum from the real data will 
be generated by dividing the observation into different epochs 
of equal duration and then cross-correlating the data cubes 
from the two epochs (Bowman et al. 2009). This approach 
preserves the persistent Hi 21 cm signal and eliminates the 
thermal noise power (which will be independent between the 
two observing epochs and, therefore, average to zero during 
the cross-correlation), leaving only the thermal uncertainty. 
Hence, the relevant noise figure for the angular power spectra 
measurement is given by: 



(_« — 
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(7) 



where A^i and N2 are simulated noise measurements from two 
different epochs (Bowman et al. 2009). 

The residual angular power spectra in the figures can be 
compared to the thermal noise uncertainty in the observations 
(Equation 7) and a fiducial 21 cm signal. We plot the ex- 
pected thermal noise uncertainty angular power spectrum of 
the MWA after 5000 hours of integration assuming a system 
temperature of Ts^.^ = 250 K, channel width of 125 KHz and 
the observing strategy described in Section 2.2. The thermal 
uncertainty spectrum is shown assuming the angular power 
spectrum has been binned in logarithmic intervals of width 
A£ = 0.1, or approximately ten bins per decade. For the ref- 
erence Hi 21 cm signal, we show the power spectrum for a 
fully neutral IGM at z = 8 (Furlanetto et al. 2006). Modeling 
the 21 cm signal using a fully neutral IGM provides a reason- 
able fiducial expectation since recent reionization simulations 
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Fig. 5. — (a) Angular power spectrum of the UVSUB residual image I"^(d,i/) made after subtraction of a foreground model with source position eiTors of 
erg = {0.01, 0.1, and 1 arcsec}. (b) Same as panel (a) but for the IMLIN residual image, I'j^i^n^id, v), that is produced after polynomial fitting and subtraction 
has been applied along each sight-line. The vertical line in panel (b) corresponds to ^ ~ 250, below which the polynomial fitting is expected to remove most of 
the structure. Both panels include the thermal noise uncertainty power spectrum assuming 5000 hours of observation with the MWA and the Hi 21 cm signal 
power spectrum for a fully neutral IGM {z ~ 8, jcm = 1; Furlanetto et al. (2006)). These angular power spectra ai'e what would be expected from the MWA if it 
integrated deep enough to directly image a typical cosmic Stromgren Sphere. 



(Lidz et al. 2008) show that the amphtude of the power spec- 
trum over the scales probed by the MWA is likely to be even 
larger than the fully neutral level when the universe if roughly 
50% ionized. It should be noted that different models predict 
different amplitudes for the Hi power spectrum. For simpli- 
fication, we have used this single realistic model to compare 
with our residual power spectrum. The specific conclusions 
regarding the scale-size dependence of where residual power 
will dominate the 2 1 cm signal will change depending on the 
reionization model. 

Figure 5(a) shows that the angular power spectrum from the 
UVSUB images are well above the thermal uncertainty power 
spectrum, as well as the model Hi 21 cm signal power spec- 
trum. In Figure 5(b), it is evident that the residual angular 



power in the IMLIN image is greatly reduced; and for two of 
our fiducial source position error levels (ae = 0.1 and 0.01 arc- 
sec), the angular power spectra intercepts the Hi signal power 
spectrum al £ < 700. This shows that the IMLIN step is very 
crucial not only for removing faint and confused continuum 
foreground sources, but also for removing residual power left 
over after subtracting the bright foreground sources. A source 
position accuracy of < 0.1 arcsec would allow the detection 
of Hi 21 cm signal at 10 < ^ < 2000 scales. 

Similar features are seen in Figure 6 for the case of calibra- 
tion errors, where the residual angular power spectrum from 
the UVSUB image only intercepts the thermal noise power 
spectrum near £ ^ 2000 and only the (Ta = 0.01 % crosses be- 
low the model Hi 21 cm signal power spectrum and the ther- 
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Fig. 6. — Same as Figure 5, but for the residuals due to calibration errors. In this case, the three residual angular power spectra are for errors of ct,, = {0.01, 0.1, 
and 1%}. 



mal uncertainty spectrum. Again, from Figure 6(b), it is evi- 
dent that the residual angular power spectrum from the IMLIN 
image is much lower, particularly below the £ = 250 threshold. 
We have not investigated in detail how scales larger than this 
threshold will be affected by the polynomial subtraction, but 
it is likely that some of the signal will be removed, as well. A 
calibration accuracy of da < 0.05% should allow the detection 
of the Hi 21 cm signal. 

3.2. ID Spherically-Averaged Power Spectrum 

The spherically-averaged three-dimensional 21 cm power 
spectrum is the primary reionization observable targeted by 
the MWA. There has been extensive research on the statisti- 
cal EoR power spectrum measurement of the brightness tem- 
perature fluctuations in low-frequency, wide-field radio ob- 
servations. Detailed formulation has been developed in the 
literature by Morales & Hewitt (2004) and Zaldarriaga et al. 



(2004). The approaches described in these efforts are inspired 
by the techniques that have been employed successfully for 
interferometric measurements of CMB anisotropics (White 
et al. 1999; Hobson & Maisinger 2002; Myers et al. 2003). 
The primary approach is to convert the full three-dimensional 
measurement cube to a one dimensional power spectrum. 

The first step is to transform our residual image cubes 
I{9, V) into y(u, rj) by performing a three dimensional Fourier 
transform denoted by the operator F({u, rf\,{9, v}). It should 
be noted here that before performing the Fourier transform, 
we have changed the units of the residual images from flux 
unit (Jy beam~^) to brightness temperature unit (mK). Hence, 
we get: 



V(vi,7^) = ¥({n,ri},{e,v})I(e,v) 



(8) 



where u,r/ = (u,v,ri). Then, we transform the measurement 
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Fig. 7. — (a) ID spherically-averaged power spectrum of the UVSUB residual image f'^d,iy) made after subtraction of a foreground model with source 
position errors of erg = {0.01, 0.1, and 1 arcsec}. (b) Same as panel (a) but for the IMLIN residual image, I"f^i^ifj(6, ^), that is produced after polynomial fitting 
and subfi'action has been applied along each sight-line. The vertical line in panel (b) con'esponds lo k = 0.03 Mpc"' , below which the polynomial fitting is 
expected to remove much of the structure. Both panels include the thermal noise uncertainty spectrum assuming 300 hours of observation with the MWA and 
binning into logarithmic spherical shells of width A.k/k = 0.5, or approximately five bins per decade. The Hi 21 cm signal power spectrum for a fully neutral 
IGM (z ~ 8, ,vhi = 1; Furlanetto et al. (2006)) is also shown. Detecting the spherically-averaged 21 cm power spectrum is the primary goal of the MWA. 



coordinates u, v into the cosmological coordinates k. 

y(k)=j(k,{u,77})y(u,77) (9) 

= J(k,{u,r,})F({u,77},{^,z.})/(0>) (10) 

where J(k, {u,r/}) denotes the Jacobian of the coordinate 
transformation from u, 77 (in units of A and //z"') to k (in units 
of cMpc~^). We have mainly followed the definition in Pee- 
bles (1993) and the formulation detailed in Morales & Hewitt 
(2004). Hence, we transformed a residual image cube (in sky 
coordinates) to a three dimensional residual visibility cube in 
the Fourier conjugate coordinates of co-moving Mpc. 

Assuming isotropy of space and ignoring redshift-space 
distortions inherent in converting our observed data cube to 



cosmological coordinates, the power spectrum can be taken 
as approximately spherically symmetric in cosmological k = 
(kjc,ky,k,) coordinates. Hence, the power spectrum can be ap- 
proximated to the square of the V(k), averaged over spherical 
shells: 

P(k)=(\V(k)\^) . (11) 

\ / \k\=k 

Thus, we obtain the one dimensional total power spectrum 
(Morales & Hewitt 2004), or the more common dimensionless 
power spectrum given by A^ = k^P{k)/(2TT^). 

While deriving the ID power spectrum, we have weighted 

the individual measurements |y(k)| by the per cell visibility 
contributions. This scheme is similar to the natural weighting 
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Fig. 8. — Same as Figure 7, but for the residuals due to calibration errors. In this case, the three residual spherically-averaged power spectra ai'e for eiTors of 
CTn = {0.01, 0.1, and 1%}. 



scheme which is applied to the raw visibilities before imaging, 
and follows the form: 



P(k): 



Eiki=.w„(k)|y(k)|^ 



E|k|=*l^«(k) 



(12) 



where Wa(k) denotes the total number of visibilities contribut- 
ing per k cell. Here, we should explicitly mention that the 
y(k) used in the above equation are the un-weighted visibili- 
ties obtained from the residual images. 

The total thermal noise power is much stronger than the ID 
spherically averaged power spectra of the Hi 21 cm signal. 
Hence, similar to the angular power spectrum case, we com- 
pare our results with the thermal noise uncertainty given by: 



P'^ik) : 



Eiku, |M(k)*A^2(k)+A^2(k)*M(k)| 



(13) 



where A^i and N2 are simulated noise measurements from two 
different epochs (Bowman et al. 2009). 

Figures 7 and 8 show the ID spherically averaged power 
spectrum from the residual images. As with the angular power 
spectrum, these figures also show theoretical Hi 21 cm power 
spectrum. However, here, instead of using a total thermal 
noise power spectrum as we did for the angular power spec- 
trum plots, we show the spherically-averaged thermal noise 
uncertainty power spectrum from 300 hours of observation 
with the MWA, as mentioned in Equation 13. The ther- 
mal uncertainty spectrum is shown assuming the spherically- 
averaged power spectrum has been binned in logarithmic 
shells of width Ak/k = 0.5, or approximately five bins per 
decade. As discussed in Lidz et al. (2008), the MWA-512 will 
be sensitive primarily to scales O.l <k < I Mpc"'. The ver- 
tical lines on Figures 7(b) and 8(b) are at A; = 0.03 Mpc"' and 
indicate the scales below which the IMLIN polynomial fit- 



10 



ting step removes significant power. These lines correspond 
to the t = 250 threshold in the angular power spectra plots. 
The higher end of the k value for the residual power spectrum 
is restricted due to the cell size in the image domain and fre- 
quency resolution of the channels in the residual image-cube. 
The maximum value of k is attained along the k. axis only, and 
hence few or no transverse (angular) modes contribute to the 
power spectrum at small scales above k > 0.6 Mpc"' in Fig- 
ures 7 and 8. The angular resolution of the MWA of ^ 4.5 ar- 
cmin (synthesized beam) corresponds to fc « 0.6 Mpc"' . 

Figure 7 shows the ID spherically-averaged power spectra 
from the residual images. These are the residual images af- 
ter the foreground subtraction in presence of source position 
errors. Figure 7(a) shows that the ID power spectra from the 
UVSUB image is well below the thermal uncertainty power 
spectrum and the Hi signal power spectrum. In Figure 7(b), 
it is evident that the residual ID spherically-averaged power 
spectra with ae = 0.01 and 0.1 arcsec from the final IMLIN 
image are below the thermal uncertainty power spectrum and 
the Hi signal power spectrum. Hence, a source position accu- 
racy of (Tg <0.l arcsec would allow the detection of Hi 21 cm 
signal with the MWA. 

Turning to the case of the calibration errors. Figure 8(a) 
shows that only the ID power spectrum from the UVSUB im- 
age for calibration error of 0.01 % is well below the thermal 
uncertainty power spectrum and the Hi signal power spec- 
trum. In Figure 8(b), it is evident that the ID spherically- 
averaged power spectra with (Ta = 0.01 and 0.1% from the 
IMLlN images are below the thermal uncertainty power spec- 
trum and the theoretical Hi signal power spectrum. Hence, 
the residual calibration accuracy of <Ta < 0.05 % would allow 
the detection of Hi 21 cm signal. 

In comparison to the angular power spectrum, we can infer 
that the ID spherically-averaged power spectra has a better 
tolerance for both the source position and residual calibration 
errors. This also reflects the fact that the angular power spec- 
trum has been produced using a single channel map of 125 
KHz, whereas the ID spherically-averaged spectrum is pro- 
duced with the total bandwidth of 32 MHz. 

3.3. Two-dimensional Power Spectrum 

In the previous section, we showed the analysis of the ID 
spherically-averaged power spectrum. However, this formu- 
lation mixes the contribution from the k i = ^ /A:? + kl and 



fell directions. It is useful, therefore, to break the averaging 
from the three dimensional k-space to the one dimensional 
A:-space into two steps since both the foregrounds and a full 
treatment of the predicted redshifted 21 cm signal that in- 
cludes redshifted-space distortions have aspherical structure 
in the Fourier domain. Following McQuinn et al. (2006), we 
first average over the transverse (angular) direction in the full 
three-dimension power spectrum to obtain P{k^, k\\). This is 
conceptually similar to averaging over the m values and keep- 
ing the i values in a CMB analysis, except we still have the 
line-of-sight dimension. Next, we obtain the 2D power spec- 
trum based on the maximum likelihood formalism following 
the same approach as used for the spherically-averaged power 
spectrum in Equations 1 1 and 12. 

Figures 9 through 1 1 illustrate the results of the simulation 
for the 2D power spectrum. In all panels of these figures, 
the color scale is held constant to facilitate comparison be- 
tween the plots. We show in Figure 9(a) a realization of the 
2D thermal noise uncertainty after 300 hours of integration 



with the MWA on our target field. Figure 9(b) shows the 2D 
Hi signal power spectrum, P(k^,k\\), of the Hi signal in units 
of mK^ Mpc"-'. Figures 10 and 1 1 show the 2D power spec- 
tra of the residual image cubes for source position errors and 
calibration errors, respectively. It should be noted that these 
plots are in different units than Figure 7 and 8. In this Sec- 
tion, we have analyzed residual images for only one of our 
fiducial error levels for each type of model corruption. For 
source position errors, we use o-g = 0. 1 arcsec and for residual 
calibration errors, we use aa =0.1%. 

Following our convention. Figure 10(a) shows that the 2D 
power spectrum from the UVSUB image, while Figure 10(b) 
shows the same for the IMLIN image. From Figure 10(a), 
we can conclude that the source position errors are more lo- 
calized towards higher k±^ values. Based on the Hi signal 
shown in Figure 9(b), we see that the signal dominates over 
the residuals at k±^ < 0.02 Mpc"'. For the IMLIN image in 
Figure 10(b), the power spectrum shows that Hi signal domi- 
nates at k±^ < 0.05 Mpc"'. The thermal uncertainty map (Fig- 
ure 9(a)) shows that thermal uncertainty dominates over the 
source position errors at lower k± values, as mentioned above. 
There are also regions with high fe|| (at k± > 0.05 Mpc"') 
where the thermal uncertainty dominates over source position 
errors. 

Figure 1 1 shows a very comparable pattern arising from the 
residual calibration errors. Hence, we find that the GSM posi- 
tion accuracy of 0.1 arcsec and calibration accuracy of 0.1% 
is sufficient to detect the Hi 21 cm signal in 2D power spec- 
trum. The advantage of the 2D spectrum over the ID spheri- 
cally averaged power spectrum is that one can even search for 
Hi signal at scales around k ^ 0.1 Mpc"' along the A:|| axes, 
which is fairly clean at lower k± values. However, in the ID 
power spectrum similar scales are dominated by the residual 
errors due to combined contribution from k± and A:|[ . 

4. CONCLUSION 

With the results from Section 3, we can address our three 
key questions for bright source subtraction residuals in the 
power spectral domain. First, we find that the level of the 
source subtraction accuracy required for power spectral de- 
tection of the 21 cm signal is roughly comparable to the accu- 
racy required for direct imaging of the Hi signal (Datta et al. 
2009). The power spectrum tolerance does suffer, however, 
compared to the direct imaging case in one regard. For direct 
imaging, it is possible to find areas in the final image map that 
are far from bright sources and have very low residual con- 
tamination. Whereas, for the power spectrum analysis, we use 
the entire image map for the calculation, mixing both the good 
and the bad areas in the image cube. Because of this differ- 
ence, the power spectrum analysis requires a more stringent 
calibration accuracy of da « 0.05% compared to da « 0.2% 
for direct imaging. This problem is most pronounced for 
the angular power spectrum analysis since the residual con- 
tamination is dominated by angular power It is mitigated 
partially in the ID spherically-averaged power spectrum by 
the inclusion of spectral information-with its lower residual 
contamination power-in the calculation, and also in the 2D 
power spectrum. However, in the case of the ID spherically- 
averaged power spectrum, it should be noted that the calibra- 
tion accuracy is determined on the basis of a ^ 300 hrs of 
integration by MWA-512. If we increase the amount of ob- 
serving time, the requirement on the accuracy of calibration 
can be lowered because our model allows the calibration error 
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Fig . 9 . — (a) Realization of the 2D thermal uncertainty power spectrum after cross-correlating simulated thermal noise maps (300 hours of total integration) 
from two different epochs (Bowman et al. 2009). (b) 2D power spectrum of the Hi 21 cm signal (Furlanetto et al. 2006) given by P(k^,k\i) = (l+2/i^+^^)P(k), 
where /i = i||/|k|. Note that the quantity plotted here and the following Figures is P{k^, ht) in units of mK^ Mpc"'. The color scale is shown in \og^QP{k_^_,k\\). 
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Fig. 10. — (a) 2D power spectrum of the UVSUB residual image / 
(7g = 0. 1 arcsec. (b) Same as panel (a) but for the IMLIN residual image, Ill^i^jj^i 
along each sight-line. The color scale is shown in \ogiQP(kj^,ku). 

to average toward zero with additional nights of observations. 

The 2-D power spectrum, in particular, addresses our sec- 
ond key question, showing clear advantages for separating the 
residual contamination from the desired signal through dis- 
tinct localization of the respective contributions in the the k± 
and A:|| plane. The results for both source position errors and 
residual calibration errors indicate that at k± < 0.05 Mpc"^ 
we are able to probe most of the A:|| scales where the Hi signal 
is dominant over the residual errors. In the ID power spec- 
trum we see dominant contribution from the residual errors 
around k ^ 0.1 Mpc"', which can be probed in the 2D power 
spectra along the k\\ axes. 

Finally, our third key question was whether we might ex- 
pect to build a template library of residual contamination er- 
rors in the power spectrum domain in order to facilitate inter- 
pretation of the final power spectrum. The results of this work 
indicate that it will indeed be possible. Further, the significant 
similarities between the 2D power spectra for the source po- 
sition error case and the gain calibration error case suggest 



, u) made after subtraction of a foreground model with source position errors of 
u), that is produced after polynomial fitting and subtraction has been applied 



that there may be common and easily identifiable properties 
of bright source residual contamination that are largely inde- 
pendent of the specific error causing the contamination. This 
will be a valuable tool for upcoming experiments. 

We have found that the IMLlN polynomial subtraction step 
is crucial not only for faint point source and diffuse fore- 
ground subtraction as studied in other works, but also for the 
success of the bright source foreground subtraction that we 
explored in this paper 

For the simulations included in this paper, we have per- 
formed the foreground subtraction of bright sources from a 
data-set of a minimum of 6 hours of observation (extrapo- 
lated to 300 and 5000 hours) in order to have the full effect 
of earth rotation synthesis. However, the MWA may per- 
form much of its bright source removal over much shorter 
timescales (^10 minutes or less) as part of its real-time cali- 
bration pipeline. The major implication for shorter time-scale 
removal of the foregrounds would be to break our assumption 
when modeling source position errors that the position errors 
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Fig. 1 1 . — Same as Figure 10, but for the residuals due to calibration errors. In this case, the residual 2D power spectrum is shown for a fiducial calibration 
error level of cto = 0. 1 %. 



are constant for the entire observation. 

We also made the assumption in this paper that each an- 
tenna's calibration errors are perfectly correlated for an entire 
6-hour observation night, but uncorrected between observ- 
ing nights. If the residual calibration error where instead per- 
fectly random between every 8-second cycle of the real-time 
calibration, then we estimate it could be possible to achieve 
the desired residual contamination noise level and detect the 
redshifted 21 cm Hi signal from reionization with a signif- 
icantly larger calibration error of da « 2.5 %. We have 
not performed our detailed simulation under this assumption, 
however, nor have we used the exact parameters that will be 
employed for the real-time calibration pipeline of the MWA. 

The results from the angular power spectrum puts more 
stringent contraints on the accuracies in source position and 
calibration. Using a larger chunk (> l25kHz)of frequency 
width might have reduced the constraints. However, we con- 
clude that if a wider bandwidth is available, it is more advan- 
tageous to perform a ID spherically averaged power spectrum 
than the angular power spectrum. This is due to the inclusion 
of the fell axes contribution to the power spectrum. 

Lastly, we would like to emphasis that similar constraints 



can also be derived for other upcoming arrays, such as LO- 
FAR and PAPER, as well as for future arrays like the Square 
Kilometer Array or a lunar array. But detailed simulations 
with the unique array specifications for each instrument would 
be required, which is beyond the scope of this paper. We ex- 
pect to build on our present analysis in future work by explor- 
ing other arrays, addressing the modified scenarios described 
above, and including additional calibration issues such as 
wide-field gain calibration of the primary beam and iono- 
sphere. 
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